clc
clear
close all
dbstop if error


dx = 0.01;
dt = 0.001;
x = -1:dx:1;  % 区域,间断在0
t = 0:dt:0.14;
index = find(x==0);

rungeKuttaStep = 3;
epsilon = 1E-10;
gamma = 1.4;
p1 = 1;
rho1 = 1;
c1 = sqrt(gamma*p1/rho1);
p2 = 0.1;
rho2 = 0.125;
c2 = sqrt(gamma*p2/rho2);
u1 = 0;
u2 = 0;
R = 8.314;

U_current = zeros(length(x),3); % 计算得到的当步的守恒变量
% 初始条件
U_current(1:index,1) = rho1;     %密度 rho
U_current(1:index,2) = rho1*u1;  %动量密度 rho*u
U_current(1:index,3) =  p1/(gamma-1) + 0.5*rho1*u1^2; %能量密度 E = p/(gamma-1) + 0.5*rho*u^2
U_current(1+index:end,1) = rho2; %密度
U_current(1+index:end,2) =  rho2*u2; %动量密度
U_current(1+index:end,3) =  p2/(gamma-1) + 0.5*rho2*u2^2; %能量密度

U1 = U_current;
U2 = U_current;
U_previous = U_current;
U_save = zeros(length(x), 3 ,length(t)); % 存储变量 质量密度，动量密度，能量密度
U_save(:,:,1) = U_current;

velocity = zeros(length(x),1);
rho = zeros(length(x),1);
pressure = zeros(length(x),1);

lamda = zeros(length(x),3);   %特征值
lamda_plus = zeros(length(x),3);
lamda_minus = zeros(length(x),3); %特征值
f_plus = zeros(length(x),3); %f+
f_minus = zeros(length(x),3); %f-
f = zeros(length(x),3); %f=(f+) + (f-)
fR_JPLUS = zeros(length(x),3); % j+1/2的右值
fL_JPLUS = zeros(length(x),3);%  j+1/2的左值
fR_JMINUS = zeros(length(x),3);
fL_JMINUS = zeros(length(x),3);


for i = 2:length(t)  % 时间循环

    for rk = 1:rungeKuttaStep
        % 首先计算每个网格点的特征值和通量
        rho = U_current(:,1);
        velocity = U_current(:,2)./rho;
        pressure = (gamma-1) * ( U_current(:,3) - 0.5*rho.*velocity.^2);
        % 特征值计算
        lamda(:,1) = velocity;
        c = sqrt(gamma*pressure./rho);
        lamda(:,2) = velocity - c;
        lamda(:,3) = velocity + c;
        lamda_plus = 0.5*(lamda + (lamda.^2 + epsilon^2).^0.5);
        lamda_minus = 0.5*(lamda - (lamda.^2 + epsilon^2).^0.5);



        f_plus(:,1) = rho/2/gamma .* ( 2*(gamma-1)*lamda_plus(:,1) + lamda_plus(:,2) + lamda_plus(:,3) );
        f_plus(:,2) = rho/2/gamma .* ( 2*(gamma-1)*lamda_plus(:,1) .* lamda(:,1) + lamda_plus(:,2) .* lamda(:,2) + lamda_plus(:,3) .* lamda(:,3) );
        f_plus(:,3) = rho/2/gamma .* ( (gamma-1).*lamda_plus(:,1) .* lamda(:,1).^2 + 0.5*lamda_plus(:,2) .* lamda(:,2).^2 + 0.5*lamda_plus(:,3) .* lamda(:,3).^2 + ...
            (3-gamma)*( lamda_plus(:,2)+lamda_plus(:,3) ).* (c.^2) /2/ (gamma-1));

        f_minus(:,1) = rho/2/gamma .* ( 2*(gamma-1)*lamda_minus(:,1) + lamda_minus(:,2) + lamda_minus(:,3) );
        f_minus(:,2) = rho/2/gamma .* ( 2*(gamma-1)*lamda_minus(:,1) .* lamda(:,1) + lamda_minus(:,2) .* lamda(:,2) + lamda_minus(:,3) .* lamda(:,3) );
        f_minus(:,3) = rho/2/gamma .* ( (gamma-1).*lamda_minus(:,1) .* lamda(:,1).^2 + 0.5*lamda_minus(:,2) .* lamda(:,2).^2 + 0.5*lamda_minus(:,3) .* lamda(:,3).^2 + ...
            (3-gamma)*( lamda_minus(:,2)+lamda_minus(:,3) ).* (c.^2) /2/ (gamma-1)); % 李老师的ppt steger-warming

        % 计算内点值
        for j = 2:length(x)-1
            if (j==2)
                %NND格式
                fR_JPLUS(j,:) = f_plus(j,:) + 0.5*minmod(f_plus(j+1,:)-f_plus(j,:) , f_plus(j,:) - f_plus(j-1,:));
                fR_JMINUS(j,:) = f_plus(j-1,:) + 0.5*minmod(f_plus(j,:)-f_plus(j-1,:) , f_plus(j-1,:) - f_plus(j-1,:));
                fL_JPLUS(j,:) = f_minus(j+1,:) - 0.5*minmod(f_minus(j+2,:)-f_minus(j+1,:) , f_minus(j+1,:) - f_minus(j,:));
                fL_JMINUS(j,:) = f_minus(j,:) - 0.5*minmod(f_minus(j+1,:)-f_minus(j,:) , f_minus(j,:) - f_minus(j-1,:));
            elseif(j==length(x)-1)

                fR_JPLUS(j,:) = f_plus(j,:) + 0.5*minmod(f_plus(j+1,:)-f_plus(j,:) , f_plus(j,:) - f_plus(j-1,:));
                fR_JMINUS(j,:) = f_plus(j-1,:) + 0.5*minmod(f_plus(j,:)-f_plus(j-1,:) , f_plus(j-1,:) - f_plus(j-2,:));
                fL_JPLUS(j,:) = f_minus(j+1,:) - 0.5*minmod(f_minus(j+1,:)-f_minus(j+1,:) , f_minus(j+1,:) - f_minus(j,:));
                fL_JMINUS(j,:) = f_minus(j,:) - 0.5*minmod(f_minus(j+1,:)-f_minus(j,:) , f_minus(j,:) - f_minus(j-1,:));

            else
                %NND格式
                fR_JPLUS(j,:) = f_plus(j,:) + 0.5*minmod(f_plus(j+1,:)-f_plus(j,:) , f_plus(j,:) - f_plus(j-1,:));
                fR_JMINUS(j,:) = f_plus(j-1,:) + 0.5*minmod(f_plus(j,:)-f_plus(j-1,:) , f_plus(j-1,:) - f_plus(j-2,:)); % j-2在第2点需要改
                fL_JPLUS(j,:) = f_minus(j+1,:) - 0.5*minmod(f_minus(j+2,:)-f_minus(j+1,:) , f_minus(j+1,:) - f_minus(j,:));% j+2在第2点要改
                fL_JMINUS(j,:) = f_minus(j,:) - 0.5*minmod(f_minus(j+1,:)-f_minus(j,:) , f_minus(j,:) - f_minus(j-1,:));
            end

            % 内点值更新
            if (rk == 1)
                U1(j,:) = U_current(j,:) - dt * (fR_JPLUS(j,:) + fL_JPLUS(j,:) - fR_JMINUS(j,:) - fL_JMINUS(j,:)) / dx;

            elseif (rk==2)
                U2(j,:) = 3/4*U_previous(j,:) + 1/4*U1(j,:) - 1/4 * dt * (fR_JPLUS(j,:) + fL_JPLUS(j,:) - fR_JMINUS(j,:) - fL_JMINUS(j,:)) / dx;

            else
                U_current(j,:) = 1/3*U_previous(j,:) + 2/3*U2(j,:) - 2/3 * dt * (fR_JPLUS(j,:) + fL_JPLUS(j,:) - fR_JMINUS(j,:) - fL_JMINUS(j,:)) / dx;
            end

        end

        if (rk == 1)
            % 边界条件
            U1(1,1) = rho1; % 给定密度
            U1(1,2) = U1(2,2); %外推
            U1(1,3) = U1(2,3);
            % 边界条件
            U1(end,1) = rho2; % 给定密度
            U1(end,2) = U1(end-1,2); %外推
            U1(end,3) = U1(end-1,3);
            U_previous = U_current;
            U_current = U1;
        end

        if (rk == 2)
            % 边界条件
            U2(1,1) = rho1; % 给定密度
            U2(1,2) = U2(2,2); %外推
            U2(1,3) = U2(2,3);
            % 边界条件
            U2(end,1) = rho2; % 给定密度
            U2(end,2) = U2(end-1,2); %外推
            U2(end,3) = U2(end-1,3);

            U_current = U2;
        end

        if (rk == 3)
 % 边界条件
            U_current(1,1) = rho1; % 给定密度
            U_current(1,2) = U_current(2,2); %外推
            U_current(1,3) = U_current(2,3);
            % 边界条件
            U_current(end,1) = rho2; % 给定密度
            U_current(end,2) = U_current(end-1,2); %外推
            U_current(end,3) = U_current(end-1,3);
            U_save(:,:,i) = U_current;
        end

    end

end
rho_final =  U_current(:,1) ;
velocity_final =  U_current(:,2)./U_current(:,1) ;
pressure_final =   (gamma-1) * ( U_current(:,3) - 0.5*rho_final.*velocity_final.^2);

rho_save = reshape(U_save(:,1,:), [] , length(t));
velocity_save =reshape( U_save(:,2,:),[] , length(t)) ./ rho_save;
pressure_save = (gamma-1) * (reshape( U_save(:,3,:) ,[] , length(t)) - 0.5 * rho_save .* velocity_save.^2 );
for j=1:length(t)

plot(x,pressure_save(:,j));
hold on
plot(x,rho_save(:,j));
plot(x,velocity_save(:,j));
legend("P","\rho","u");
drawnow;
pause(0.1);
hold off
end